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Abstract. This article presents a general and novel approach to the automation of goal-oriented 
error control in the solution of nonlinear stationary finite element variational problems. The approach 
is based on automated linearization to obtain the linearized dual problem, automated derivation and 
evaluation of a posteriori error estimates, and automated adaptive mesh refinement to control the 
error in a given goal functional to within a given tolerance. Numerical examples representing a 
variety of different discretizations of linear and nonlinear partial differential equations are presented, 
including Poisson's equation, a mixed formulation of linear elasticity, and the incompressible Navier- 
Stokes equations. 
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1. Introduction. For any numerical method, it is of critical importance that 
the accuracy of computed solutions may be assessed. For the numerical solution of 
differential equations, the accuracy is typically assessed manually by computing a 
sequence of solutions using successive refinement until it is judged that the solution 
has "converged". This approach is unreliable as well as time-consuming. It may also 
be impossible, since computing resources may be exhausted long before convergence 
has been reached. 

For finite element discretizations, classic a posteriori error analysis provides a 
framework for controlling the approximation error measured in some Sobolev norm, 
cf. [l]. Over the last two decades, goal-oriented error control has been developed as 
an extension of the classic a posteriori analysis [9l |12] . The problem of goal-oriented 
error control for stationary variational problems can be posed as follows. Consider 
the following canonical variational problem: find u such that 

F{u]v)={) MveV, (1.1) 

where FiFxt^^Risa semilinear form (linear in v) on a pair of trial and 
test spaces (V, V"). A goal-oriented adaptive algorithm seeks to find an approximate 
solution Uh ^ u oi such that 

r]=\M{u)-M{uh)\ <e, 

where 7W : F ^ M is a given goal functional, e > is a given tolerance, and 77 is here 
defined as the error in the given goal functional. In other words, goal-oriented error 
control allows the construction of an adaptive algorithm that targets a simulation to 
the efficient computation of a specific quantity of interest. 

The framework developed in [9l |T2] provides a general method for deriving an 
a posteriori estimate of the error and adaptive refinement indicators, based on the 
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solution of an auxiliary linearized adjoint (dual) problem. This framework is directly 
applicable to a large class of finite element variational problems. However, a certain 
level of expertise is required to derive the error estimate for a particular problem and 
to implement the corresponding adaptive solver. In particular, the derivation of the 
dual problem involves the linearization of a possibly complicated nonlinear problem. 
Furthermore, both the derivation and evaluation of the a posteriori error estimate 
remain nontrivial (at least in practice). Moreover, each derivation must be carried 
out on a problem- by-problem basis. As a result, goal-oriented error control remains 
a tool for experts and usually requires a substantial effort to implement. 

In this work, we seek to automate goal-oriented error control. We present a fully 
automatic approach for computation of error estimates and adaptive refinement; that 
is, without any manual analysis, preparation, or intervention. To our knowledge, no 
such automation has previously been presented, and in particular not realized. The 
strategy presented here requires a minimal amount of input and expert knowledge; 
the only input required by our adaptive algorithm is the semilinear form F, the func- 
tional A^, and the tolerance e. Based on the given input, the adaptive algorithm 
automatically generates the dual problem, the a posteriori error estimate, and at- 
tempts to compute an approximate solution that meets the given tolerance for 
the given functional. In particular, problem-tuned error estimates and indicators are 
generated without any manual derivations. This has the potential of rendering state- 
of-the-art goal-oriented error control fully accessible to non-experts at no additional 
implement ational cost. 

We emphasize that although one may, in principle, manually carry out the nec- 
essary analysis and implementation in any particular case, automation plays an im- 
portant role since it (i) makes expert knowledge accessible to non-experts, (ii) speeds 
up development cycles, and (iii) enables more complex problems to be tackled which 
would otherwise require considerable effort to analyze and implement. A similar and 
successful automation effort as part of the FEniCS Project [271 Ull IlQl HI] has led to 
the development of methodology and tools that automate the discretization of a large 
class of partial differential equations by the finite element method. 

In this paper, we limit the discussion to stationary variational problems. The 
error estimates and indicators generated by the automated algorithm can be viewed 
as a version of the dual- weighted-residual estimates of [9 . We emphasize that the 
main target of this paper is to present an automation of a method for goal-oriented 
error control; in particular, it is not our intention to improve the method theoretically 
nor in detail examine the cases where the method is known to work poorly. Also, the 
question of automated error control for time-dependent problems will be considered 
in later works. 

The remainder of this introduction describes the organization of the paper. The 
first two sections. Sections [2] and [Sj establish the abstract problem setting by defining 
notation and summarizing the well-established goal-oriented error estimation frame- 
work for linear variational problems. The linear setting is discussed first for the sake 
of clarity and brevity: the extension to the nonlinear case is summarized in Section [6j 
The primary novel contributions of this paper are contained in Sections |4j [5) and [7| 
the main result of this paper is an automated strategy for the computation of error 
estimates and indicators. This strategy relies on two key components: first, a proce- 
dure, applicable to a general class of stationary variational problems, for the derivation 
of a strong residual representation from a weak residual representation (Section [4|. 
Second, the evaluation of the error estimates relies on a dual approximation: an af- 



AUTOMATED GOAL-ORIENTED ERROR CONTROL I 



3 



fordable strategy for obtaining an improved dual approximation by extrapolation for 
a general class of finite element spaces is described in Section [5] 

Combining the goal-oriented error estimation framework with the automated tech- 
niques introduced in Section |4] and |5j yields the automated adaptive algorithm de- 
scribed in Section [71 A realization of this automated algorithm has been implemented 
as part of the FEniCS Project, and is, in particular, available through both the 
Python and C++ versions of the DOLFIN library. A simple example of its use and 
some aspects of the implementation are also discussed in this section. In Section [Sj we 
apply the presented framework to three examples of varying complexity: the Poisson 
equation, a three-field mixed formulation for the linear elasticity equations, and the 
stationary Navier-Stokes equations. Finally, we conclude and discuss further work in 
Section H 

2. Notation. Throughout this paper, ^7 C denotes an open, bounded domain 
with boundary dQ. We will generally assume that Q is polyhedral such that it can be 
exactly represented by an admissible, simplicial tessellation Th- The boundary will 
typically be the union of two disjoint parts, denoted dClo and 917 at- 

In general, the notation V{X]Y) is used to denote the space of fields X ^ Y 
with regularity properties specified by V. If Y = R, this argument is omitted. For 
M^); that is, the space of (i- vector fields on C 1] in which each component 
is square integrable, the inner product reads (•, and the norm is denoted || • Wk- 
If K = the subscript is omitted. For m = 1,2,..., H'^{Q) denotes the space 
of square integrable functions with m square integrable distributional derivatives. 
Also, H^Y = {'^ ^ H^{Q) : u\y = g}. Similarly, i^(div,r^) denotes the space of square 
integrable vector fields with square integrable divergence. Note that both the gradient 
of a vector field and the divergence of a matrix field are applied row- wise. 

A form a :Wi X ■ ■ ■ X Wn x Vp x ■ ■ ■ x Vi ^ written a{wi, . . . , Wn; Vp,. . . , vi), 
is (possibly) nonlinear in all arguments preceding the semi-colon, but linear in all 
arguments following the semi-colon. 

3. A framework for goal-oriented error control. In this section, we present 
a general framework for goal-oriented error control for conforming finite element dis- 
cretizations of stationary variational problems. The framework is a summary of the 
paradigm developed in [9l [12] . For clarity, we restrict our attention to linear vari- 
ational problems and linear goal functionals. Extensions to nonlinear problems and 
nonlinear goal functionals are made in Section |6] 

Let V and V be Hilbert spaces of functions or fields defined on a domain 1] C 
for (i = 1,2,3. In this section, we consider the following linear variational problem: 
find u such that 

a{u,v)=L{v) MveV. (3.1) 

We assume that a : V xV ^ R is continuous, bilinear form, and that L : V ^ R 
is a continuous, linear form. We shall further assume that the problem is well-posed; 
that is, there exists a unique solution u that depends continuously on any given data. 



The variational problem defined by (3.1) will be referred to as the primal problem 
and u will be referred to as the primal solution. 

Let Th be an admissible simplicial tessellation of ^ (to be determined) and assume 
that Vh C V and Vh C V are finite element spaces defined relative to Th- The finite 



element approximation of (3.1) then reads: find G Vh such that 



a{uh,v) = L{v) yveVh- (3.2) 
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We assume that the spaces Vh and Vh satisfy an appropriate discrete inf-sup condition 



such that a unique discrete solution exists. The problem (3.2) will be referred to as 



the discrete primal problem and Uh the discrete primal solution. 

We are interested in estimating the magnitude of the error in a given goal func- 
tional M. : F ^ M. Moreover, for a given tolerance e > 0, we aim to find (Vh^Vh) such 
that the corresponding finite element approximation as defined by (3.2), satisfies 



T] = \M{u) - M{uh)\ <e. 



(3.3) 
effi- 



In addition, we would like to compute the value of the goal functional A4{uh 
ciently, ideally using a minimal amount of work. 

In order to estimate the magnitude of the error r^, we first define the (weak) 
residual relative to the approximation u^, 



r{v) = L{y) - a{uh,v). 



(3.4) 



Some remarks are in order. First, r is a bounded, linear functional by the continuity 
and linearity of a and L. Second, as a consequence of the Galerkin orthogonality 
implied by V/^ C V", the residual vanishes on Vh- In other words. 



r{v)=0 yveVh. 
Next, we define the (weak) dual problem: find z G such that 
a*(z,v) =M{v) yv e T>*, 



(3.5) 



(3.6) 



where (V*, V"*) is the pair of dual trial and test spaces, and a* denotes the adjoint of a; 
that is, a*('U, w) = a{w^ v). We shall assume that the dual problem (3.6) is well-posed. 



and that there thus exists a dual solution z solving (3.6) with continuous dependence 



on the input data. Moreover, we assume that the dual trial and test spaces are chosen 
such that u - Uh e V* and z e V. This holds if V'' = Vq = {v - w : v,w e V} and 
V* = V. 



Combining (3.6), (3.4), and (3.1), we find that 



M{u) - M{uh) = a*(z, u-Uh) = a{u - uh, z) = L{z) - a{uh, z) = r{z). 

The error M{u) — M{uh) is thus equal to the (weak) residual r evaluated at the 
dual solution z. By the Galerkin orthogonality (3.5), we obtain the following error 
representation: 



M(u) - M{uh) = r{z) = r{z - tthz). 



(3.7) 



Here, iThZ G Vh is an arbitrary test space field, typically an interpolant of the dual 
solution. 

An identical error representation is obtained for nonlinear variational problems 
and nonlinear goal functionals with a suitable definition of the dual problem. We re- 
turn to this issue in Section [6] It follows that if one can compute (or approximate) the 
solution of the dual problem, one may estimate the size of the error by a direct eval- 
uation of the residual. However, some concerns remain that require special attention. 



First, the error representation (3.7) is not directly useful as an error indicator. The 



derivation of an a posteriori error estimate and corresponding error indicators from 
the error representation has traditionally required manual analysis, typically involving 
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some form of integration by parts and a redistribution of boundary terms (fluxes) over 



cell facets. Second, for (3.7) to give a useful estimate of the size of the error, care must 



be taken when solving the dual problem (3.6). In particular, the error representation 



evaluates to zero if the dual solution is approximated in 14 . Finally, the derivation of 
the dual problem may involve the differentiation of a nonlinear variational form. We 
discuss how each of these issues can be automated in the subsequent sections. 

4. Automated derivation of error estimates and error indicators. This 
section presents a novel approach to the derivation of an a posteriori error estimate 
and corresponding error indicators for a general class of stationary variational prob- 



lems. The starting point is the general abstract error representation (3.7). The 
proposed generic representation for the estimate and indicators is motivated and in- 
troduced in Section |4.1[ followed by a strategy for computing this representation in 



Sections |4.2f|4.3[ The key idea is the computation of a problem-tuned strong residual 
representation using only ingredients from the general abstract form; this concept 
constitutes one of the pillars for the complete automated strategy. We emphasize 
that the resulting error estimate coincides with classical duality-based error estimates 
that may be derived manually (using integration by parts) for standard problems 
such as the Poisson problem. However, we present here an approach that allows these 
error estimates to be generated and evaluated automatically. We also note that the 
automatically derived error indicators differ from the standard duality-based error 
estimates resulting from an integration by parts followed by one or more inequalities. 

4.1. A generic residual representation. For motivational purposes, we start 
by considering the standard derivation of a dual-weighted residual error estimator for 
Poisson's equation: — Aii = / and its corresponding variational problem defined by 
a{u^v) = (gradix, grad'u) and L{v) = {f^v) on V = V = Hq{Q). By integrating the 
weak residual by parts cell- wise, one obtains 

r{z) = L{z) - a{uh,z) = {f,z) - (gradix/^, grad2:) 

= ^{f,z)T - {graduh^grad z)t ^ (/ + Auh, z)t + {-dnUh, z)ot 
TeT TeT 

= ^ (/ + Aii/,, z)t + {[-dnUh],z)dT^ 

TeT 

where [—dnUh] denotes an appropriate redistribution of the flux over cell facets. 
Several choices are possible, see for example [1, Chap. 6], but we here make the 
simplest possible choice and distribute the flux equally. In particular, we deflne 
[Snii/i]!^ = ^{grdiduhlr ' ^ + gradii/i|T' • n') over all internal facets S shared by two 
cells T and T\ and [dnUh]\s = dnUh\s on external facets (facets on the boundary of 
Q). Hence, one may estimate the error by 

\M{u)-M{uh)\ < ^7?T, (4.1) 

TGT 

where the error indicator tjt is given by 

VT = \{f + AUh.Z- 7lhZ)T + {[-dnUh],Z- 7lhZ)dT\' (4.2) 

We note that although one may in principle use rjT = | (/, — (grad u^^ grad z)t\ 
as an error indicator (without integrating by parts and redistributing the normal 
derivative), that indicator is much less efficient than the error indicator defined 
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in (4.2). Both indicators will sum up to the same value (if taken with signs), but 



only as a result of cancellation. The error indicator (4.2) is generally smaller in mag- 
nitude, scales better with mesh refinement, and gives a sharper error bound when 
summed without signs. See for an extended discussion. 

Estimates similar to ( |4.1[ ) have been derived by hand (originally for use with 
norm-based error indicators) for a variety of equations. A non-exhaustive list of 
examples (including purely norm-based indicators) includes standard finite element 
discretizations of the Poisson equation [3], various mixed formulations for the Stokes 
equations and stationary Navier-Stokes equations ^37|, i7(div)-based discretizations 
of the mixed Poisson and mixed elasticity equations [H|30], and i^(curl)-based dis- 
cretizations for problems in electromagnetics [7J. Duality-based goal-oriented error 
estimates have been derived for a number of applications, including ordinary differ- 
ential equations [13], plasticity [32], hyperbolic systems [22], reactive compressible 
flow [36] , systems of nonlinear reaction-diffusion equations [M] |35] , eigenvalue prob- 
lems pT, wave propagation [4j, radiative transfer [33 , nonlinear elasticity ^25j, the 
incompressible Navier-Stokes equations [H |T8] , variational multiscale problems [24] , 
and multiphysics problems [23]. 

These estimates share a common factor, namely that the error is expressed as a 
sum of contributions from the cells and the facets of the mesh. Moreover, each of these 
estimates has been derived manually for the specific problem at hand. Here, we aim 
to demonstrate that for a large class of variational problems, one may automatically 
compute an equivalent residual representation. The representation takes the following 
generic form: 

r{v) = Yl (^T:V)t + {RdT:V)QT = {Rt.v)t + [{RdT.v)QT], (4.3) 
TeTh TeTh 

where 

[{RdT,v)dT]= Y \{i^dT,v\T)s ^ {RdT',v\T')s) ^ ^ {RdT,v)s- 
SCdTnQ SCdTndQ 

It follows that one may use as error indicators 

Vt = \ {Rt^z -7rhz)T ^ [{Rqt.z -7rhz)QT]\- (4.4) 

In these expressions, Rt denotes a residual contribution evaluated over the domain of a 
cell T, whereas Rqt denotes a residual contribution evaluated over a cell boundary dT. 

4.2. Automatic computation of the residual representation. We shall 
focus our attention on a class of residuals r satisfying the following assumptions: 
Al (Global decomposition) The residual is a sum of local contributions: 

r = ^ TT. 

TeTh 

A2 (Local decomposition) Each local residual offers a local decomposition: 

rriv) = {Rt, v)t + {RdT, v)dT V v G V\t- (4.5) 



We note that Al is satisfied if the bilinear and linear forms a and L are expressed 
as integrals over the cells and facets of the tessellation Th- We also note that A2 is 
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Fig. 4.1. The bubble function br- 



satisfied if the variational problem (3.1) has been derived by testing a partial differ- 



ential equation against a test function and (possibly) integrating by parts to move 
derivatives onto the test function. 

For the sake of a simplified analysis, we also make the following assumption: 
A3 (Polynomial representation) The residual contributions (or, in the case of vector 
or tensor fields, each scalar component of these) are piecewise polynomial: 

We discuss the implications of this assumption below, but note that one of the nu- 
merical examples presented does not satisfy this assumption. 

In the following, we shall show that a residual representation may be automatically 
computed if assumptions A1-A2 are satisfied. More precisely, if assumptions A1-A3 
are satisfied, we shall show that one may automatically compute the exact residual rep- 



resentation (4.3) for a given variational problem (3.1). In particular, one may directly 



compute the cell and facet residuals Rt and Rqt by solving a set of local problems on 
each cell T of the tessellation Th- If A3 fails; that is, if only assumptions A1-A2 are 
satisfied, the automated procedure computes weighted L^-projections of the residual 
decomposition terms and hence an approximate residual representation. 

To compute the cell residual Rt^ let be a basis for V^{T) and let 6t 

denote the bubble function on T. We recall that for a simplex T C M^, the bubble 
function hr is defined by 



where A^. is the barycentric coordinate function on T associated with vertex Xi (the 
iih linear Lagrange nodal basis function on T). Note that 6t vanishes on the boundary 



of T. See Figure [471] for an illustration. Testing the local residual tt against 6^^*, 
we obtain the following local problem for the cell residual Rt'- find Rt G V^{T) such 
that 



{Rt, bT^i)T = rT{hT(l)i), i = 1, . . . m. 



(4.6) 



To obtain a local problem for the facet residual RdT, we define for each facet S 
on T the cone function /3j by 



T 



(4.7) 
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Fig. 4.2. The cone function (3^ . 



where /J is a suitably defined inde x set such that 



For an ihustration, see Figure 



4.2 



on ah facets / of T but S. 



Clearly, /3^\s = bs- Next, let {(j)i}2=i be a basis 



for V^{T). Testing the local residual vt against we obtain the following local 

problem for each facet residual: find RdT\s ^ V^{S) such that 



Vi G /J. 



(4i 



We prove below that by assumptions A1-A3, the local problems (4.6) and ( |4.8[ ) 
uniquely define the cell and facet residuals Rt and Rqt of the residual represen- 
tation 



One may thus compute the residual representation (4.3) by solving a set of local 
problems on each cell T. First, one local problem for the cell residual Rt^ and 
then local problems for the facet residual Rqt restricted to each facet of T. If the 
test space is vector- valued, the local problems are solved for each scalar component. 



We emphasize that the computation of the residual representation (4.3) and thus 



the error indicator (4.4) may be computed automatically given only the variational 



problem (3.1 ) in terms of the pair of bilinear and linear forms a and L. In particular. 



the derivation of the error indicators does not involve any manual analysis. 

We remark that the use of bubble and cone functions to localize the weak residual 
is a standard technique for proving reliability and efficiency for norm-based error 
estimates (see e.g. [38]). The crucial observation here is that this technique can also 
be used to automatically generate the residual decomposition given only the weak 



residual. We also remark that the local problems (4.6) and (4.8) are different from the 



local problems that were introduced in [6 to represent the cell and facet residuals Rt 
and Rqt as a single residual. 

4.3. Solvability of the local problems. To prove that the local problems (|4.6|) 



and (4.8) uniquely determine the cell and facet residuals, we recall the following result 
regarding bubble- weighted L^-norms. For a proof, we refer to [TJ Theorems 2.2, 2.4]. 

Lemma 4.1. Let T he a d-simplex and let br denote the bubble function on T. 
There exist positive constants c and C , independent of T, such that 



cM\l<{hT4>,' 



)t<cm\1 



(4.9) 



for all4>erP{T). 

We may now prove the following theorem. 
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Theorem 4.1. // assumptions A1-A3 hold, then the cell and facet residuals of 



and (4.8). 



the residual representation (4.3) are uniquely determined by the local problems (4.6) 



Proof Consider first the cell residual Rt- Take v = bTcj^i in (4.5) for i = 1, . . . , m. 
Since v vanishes on the cell boundary 9T, we obtain (4.6). By assumption, Rt G 



V^{T) and is thus a solution of the local problem (4.6). It follows from Lemma 4.1 
that it is the unique solution. We similarly see that the facet residual Rqt is a solution 
of the local problem (4.8) and uniqueness follows again from Lemma 4.1 □ 



In the cases where A3 fails, such as if the variational problem contains non- 



polynomial data, the local problems (|4.6|) and (4.8) uniquely determine the projections 



of Rt and RotIs onto V^{T) and V^{S) respectively. The accuracy of the approxi- 
mation may then be controlled by the polynomial degrees p and q. In the numerical 
examples presented in Section [Sj we lei p = q be determined by the polynomial degree 
of the finite element space. We have not observed any significant errors introduced 
by this approximation in our numerical experiments. 

Looking back at the special case of Poisson's equation (4.2), the cell and facet 
residuals derived by hand are given by Rt = / + and Rqt = —dnUh^ respec- 
tively. We emphasize that, by what is just shown, this indeed coincides with the 



representation defined by (4.6) and (4.8) if / is polynomial. 



5. Approximating the dual solution. In order to evaluate the error repre- 
sentation (3.7) and to compute the error indicators (4.4), one must compute, or in 



practice approximate, the solution z of the dual problem (3.6). The natural discretiza- 
tion of (3.6) reads: find = Vh such that 



a*{zh,v)=M{v) ^veV:=Vh,o- 



(5.1) 



However, since the residual r vanishes on Vh^ Zh is, for the purpose of error estimation, 
highly unsuitable as an approximation of the dual solution. 

An immediate alternative is to solve the dual problem using a higher order 
method. If the dual solution is sufficiently regular, a higher order method would 
be expected to give a more accurate dual approximation. It is observed in practice 
that a more accurate dual approximation gives a better error estimate [9 , although 
complete reliability cannot be guaranteed [31 . Other alternatives include approxima- 
tion by hierarchic techniques [l|i5j or approximating the dual problem on a different 
mesh. In this work, we suggest a new alternative based on solving (5.1) using the 



same mesh and polynomial order as the primal problem and then extrapolating the 
computed solution Zh to a higher order function space. This strategy can be compared 
to the higher order interpolation procedure presented in [9 for regular quadrilater- 
al/hexahedral meshes. The strategy presented here extends that of [9 however, as it 
can be applied to almost arbitrary (admissible) simplicial tessellations. 

To define the extrapolation procedure, let Vh be a finite element space on a tes- 
sellation Th and let Wh D 14 be a higher order finite element space on the same 
tessellation Th- Furthermore, let {(/>J}j^^i be a local basis for Wh on T and let 
{^iljLi be the corresponding global basis. For Vh G 14, we define the extrapola- 
tion operator E : Vh ^ Wh as described in Algorithm [l] This algorithm computes 
the extrapolation by fitting local polynomials to the finite element function Vh on 
local patches. This yields a global multi- valued function which is then averaged to 



obtain the extrapolation Evh- We illustrate the extrapolation algorithm in Figure 5.1 
for a one-dimensional case. 
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Algorithm 1 Extrapolation 

1. (Lifting) For each cell T eTh' 

(a) Define a patch of cells uot ^ T oi sufficient size and let be the 
collection of degrees of freedom for Vh on the patch. The size of the 
patch ujt should be such that the number of degrees of freedom m is 
greater than or equal to the local dimension n of 1V^|t- 

(b) Let be a smooth extension of {^Jj^^^i to the patch ujt- 

(c) Define Aij = £^(0^^) and bi = ii{vh) for z = 1, . . . , m, j = 1, . . . , n. 

(d) Compute the least-squares approximation of the (overdetermined) 
m X n system A^t = b. 

2. (Smoothing) 

(a) For each global degree of freedom j, let Xj be the set of corresponding 
local expansion coefficients determined on each cell T by the local vector 
^T- Define = '\Y~\'^xeXj ^' ^^^^ ^^^^ ^ ^ degrees of 
freedom that are shared between cells. 

(b) Define Evh = EjLi O^i- 



Algorithm [T] may be used to compute a higher order approximation of the dual 
solution z as follows. First, we compute an approximation G Vh of the dual solution 
by solving ( |5.1| ). We then compute the extrapolation Ezh G Wh where Wh is the finite 
element space on Th obtained by increasing the polynomial degree by one. We then 
estimate the error by 

V = \M{u) - M{uh)\ = \r{z)\ ^ \r{Ezh)\ = Vh- 

6. Extensions to nonlinear problems and goal functionals. We now turn 
to consider nonlinear variational problems and goal functionals. We consider the 
following general nonlinear variational problem: find u G V such that 

F{u',v)=0 yveV. (6.1) 

For a given nonlinear goal functional : F ^ R, we define the following dual 
problem: find z G V* such that 

F*(z,v) = Af(v) yv e T>*, (6.2) 

where, as before, = Vq and V* = V . The bilinear form F' is an appropriate 
average of the Frechet derivative F'{u] 5u^ v) = ^^^^^ Su of F, 

F^{-r)= f F\su^{l-s)uhrr)ds. (6.3) 

We note that by the chain rule, we have F^{u — Uh-,-) — F{u; •) — F{uh] •)• The linear 
functional M' is defined similarly. Note that ( |6.2| ) reduces to ( |3.6| ) in the linear case 
where F{u] v) = a{u^ v) — L{v). 

The following error representation now follows directly from the definition of the 
dual problem: 

M{u) - M(uh) = M'(u - Uh) = F'*{z, u-Uh) = F'(u - Uh, z) 
= F{u] z) - F{uh] z) = -F{uh; z) = r{z). 
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Fig. 5.1. Extrapolation of a continuous piecewise linear function Vh to a continuous piecewise 
quadratic function Ev^. The extrapolation is computed hy first fitting a quadratic polynomial on each 
patch. In one dimension, each patch is a set of three intervals and each local quadratic polynomial 
is computed by solving an overdetermined 4x3 linear system. The continuous piecewise quadratic 
extrapolation Evh is then computed by averaging at the end-points of each interval. 



We thus recover the error representation (3.7). 



In practice, the exact solution u is not known and must be approximated by 
the approximate solution Uh] that is, the linear operator F' is approximated by the 
derivative of F evaluated ai u = u^. The resulting linearization error may for the sake 
of simplicity be neglected, as we shall in this exposition, but doing so may reduce the 
accuracy (and reliability) of the computed error estimates. For a further discussion on 
the issue of linearization errors in the definition of the dual problem, we refer to [9] . 

It follows that the techniques described in Section [4] and [5] directly apply to the 
residual r and the dual approximation Zh also for the nonlinear case. 



7. A complete algorithm for automated goal-oriented error control. 

Based on the above discussion, we may now phrase the complete algorithm for auto- 
mated adaptive goal-oriented error control in Algorithm [2] 
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Algorithm 2 Adaptive algorithm 



Let F : V X V ^ M he di given semilinear form, let A1 : F ^ M be a given goal 
functional, and let e > be a given tolerance. 

1. Select an initial tessellation Th of the domain Q and construct the correspond- 
ing trial and test spaces Vh C V and V/i C F (for a given fixed finite element 
family and degree). 



Compute the finite element solution Uh G Vh of the primal problem (6.1) 
satisfying F{uh]v) =0 for all v G V^. 



3. Compute the finite element solution G of the dual problem (6.2) satis- 
fying v) = M'{v) for ah v e . 

4. Extrapolate ^ Ez^ using Algorithm [l] 

5. Evaluate the error estimate rjh = \F{uh'-, Ezh)\- 

6. If l^^^l < e, accept the solution Uh and break. (Stopping criterion) 

7. Compute the cell and facet residuals Rt and Rqt of the residual representa- 
tion ( |4.3| ) by solving the local problems ( |4.6[ ) and ( |4.8[ ). 

8. Compute the error indicators 

VT = \ {RT,EZh - 7lhEZh)T + [{RdT.EZh - 7rhEZh)dT]\' 

9. Sort the error indicators in order of decreasing size and mark the first M 
cells for refinement where M is the smallest number such that VTi ^ 
a J^TeTh some choice of a G (0, 1]. (Dorfier marking [11]) 

10. Refine all cells marked for refinement (and propagate refinement to avoid 
hanging nodes). 

11. Go back to step 2. 

Algorithm [2] has been implemented within the FEniCS project [27l [26] [29], a 
collaborative project for the development of concepts and software for automated 
solution of differential equations. The implementation is freely available, and dis- 
tributed as part of DOLFIN (version 0.9.11 and onwards). We discuss some of the 
features of the implementation here and provide a simple use case. More details of 
the implementation will be discussed in future work [34] , 

For the specification of variational problems, the Python interface of DOLFIN 
accepts as input variational forms expressed in the form language UFL [2 . Forms 
expressed in the UFL language are automatically passed to the FEniCS form com- 
piler FFC [20] [111 [28] which generates efficient C++ code for finite element assembly 
of the corresponding discrete operators. For a detailed discussion, see [29^. Stationary 
discrete variational problems can be solved in DOLFIN by calling the solve function 
accepting as input a variational problem specified by a variational equation expressed 
by two variational forms (defining the left- and right-hand sides), the solution func- 
tion u and any boundary conditions bcs. Our implementation adds the possibility 
of solving such problems adaptively with goal-oriented error control by adding a goal 
functional M and an error tolerance, say le-6: 



solve (a 


== L, u. 


bcs , 


tol = l 


e-6 , 


M = M) 


# Linear case 


solve (F 


== 0, u. 


bcs , 


tol = l 


e-6 , 


M = M) 


# Nonlinear case 



A simple complete example is listed in Figure 7.1 A number of optional parameters 
may be specified to control the behavior of the adaptive algorithm, including the 
marking strategy and the refinement fraction. The default marking strategy is Dorfier 
marking [11 with a refinement fraction of a = 0.5. 

Internally, the adaptive algorithm relies on the capabilities of the form Ian- 
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from dolfin import * 

mesh = Unit Square (4 , 4) 

V = Funct ionSpace (mesh , "CG", 1) 
u = Function(V) 

V = TestFunct ion (V) 
f = Constant ( 1 . 0) 

F = inner((l + u* * 2 ) * grad (u) , grad(v))*dx - f*v*dx 
be = DirichletBC (V , 0.0, "near(x[0], 0.0)") 

M = u*dx 

solve(F == 0, u, be, tol=l.e-3, M=M) 



Fig. 7.1. Complete code for the automated adaptive solution of a nonlinear Poisson-like problem 
on the unit square with f = 1.0, homogeneous Dirichlet boundary conditions on the left boundary and 
homogeneous Neumann conditions on the remaining boundary, with goal functional A4 = f^udx. 



guage UFL for generating the dual problem, the local problems for the cell and facet 
residuals, and the computation of error indicators. As an illustration, we show here 
the code for generating the bilinear form a* = F^* of the dual problem (6.2): 



a_star = ad j oint ( der ivat ive (F , u)) 



8. Numerical examples. In this section, we aim to investigate the performance 
of the automated algorithm. Since the theoretical properties of the proposed extrapo- 
lation procedure are largely unknown, the investigation here focuses on the quality of 
the error estimate and the sum of the error indicators on adaptively refined meshes. 
The total computational efficiency of the automated adaptive algorithm will be inves- 
tigated in later works [34] . 

We present three numerical examples from three different application areas, aim- 
ing to illustrate different characteristics and varying levels of complexity. We begin 
by considering a basic example: a standard discretization of the Poisson equation, 
and evaluate the quality of the error estimates; the results show that the algorithm 
gives error indicators close to the optimal value of one. The second example is a 
discretization of a weakly symmetric formulation for linear elasticity. This discretiza- 
tion involves a nontrivial finite element space, namely, a mixed finite element space 
consisting of multiple Brezzi-Douglas-Marini elements, and multiple discontinuous 
and continuous elements. As far as the authors are aware, this is the ffist demon- 
stration of goal-oriented error control for the discretization presented. The results 
show that the algorithm produces error estimates of optimal quality also for this far 
more complicated case. Finally, we consider a nonlinear, nonsmooth example of wide- 
spread use: a mixed discretization of the incompressible Navier-Stokes equations and 
evaluate both the quality of the error estimates and the performance of the adaptive 
algorithm. 

8.1. The Poisson equation. We begin by considering the Poisson equation: 

-Au = / in (8.1a) 
u = OondnD, (8.1b) 
dnU = g on dClN- (8.1c) 
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The standard variational formulation of (8.1) fits the framework of Section [s] with 
V = V = Hlg^jn)and 

a(ii, 'u) = (gradii, gradi'), (8.2a) 
L{v) = {f,v)^{g,v)9nr.- (8.2b) 



We consider the discretization of ( |8.2[ ) using the space of continuous piecewise linear 
polynomials that satisfy the essential boundary condition for Vh = Vh- 
As a test case, we consider a three-dimensional L-shaped domain, 

n = {{-i,i) X (-1,1) \ (-1,0) X (-1,0)) X (-1,0), 

with Dirichlet boundary OVLd = {(^7^7^) • x = 1 01 y = 1} and Neumann boundary 
OVLn = dft \ dVtD' Let /(x, z) = —2{x — 1) and let = G • n with G(x, y^ z) = 
{{y — 1)^, 2(x — l){y — 1), 0). The exact solution is then given by 

u{x,y,z) = {x-l){y-l)\ (8.3) 

As a goal functional, we take the average value of the solution on the left boundary 
r = {{x^y^z) : X = —1}; that is, 

M{u) = J uds. (8.4) 
It follows that the exact value of the goal functional is M{u) = —2/3. 



Figure 8.1 shows errors 77, error estimates 77/^, the sum of the error indicators 
Vt^ and efficiency indices r]h/v and tjt/t] for a series of adaptively (and auto- 
matically) refined meshes. We first note that the error estimate r]h is very close to the 
error 77. On the coarsest mesh, the efficiency index is 77/1/77 ~ 0.89 and as the mesh 
is refined, the efficiency index quickly approaches r]h/r] ^ 1. We further note that 
the sum of the error indicators tends to overestimate the error, but only by a small 
constant factor. This demonstrates that the automatically computed error indicators 
are good indicators for refinement. We emphasize that since the error indicators are 
not used as a stopping criterion for the adaptive refinement, it is not important that 
they sum up to the error. 

8.2. Weakly symmetric linear elasticity. As a more challenging test prob- 
lem, we consider a three-field formulation for linear isotropic elasticity enforcing the 
symmetry of the stress tensor weakly. This gives rise to a mixed formulation that 
involves i^(div)- and L^-conforming spaces. For a domain C M^, the unknowns 
are the stress tensor a G i^(div, 1^; R^^^), the displacement u G L^(r^;R^), and the 
rotation 7 G L'^{Q). The bilinear and linear forms read 

«((c^,^,7), (^,^,^)) = {Act.t) + (diver, + (ii,divr) + (a, 7?) + (7,r), (8.5a) 
L((r, V, T])) = {g, v) + {uq, t • n)d^i- (8.5b) 

Here, ^ is a given body force, uq is a prescribed boundary displacement field, and 
A is the compliance tensor. For isotropic, homogeneous elastic materials with shear 
modulus ji and stiffness A, the action of A reduces to 

^^=^i^- . (8.6) 

2/i V 2(/i + A) ^ V 
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(b) Efficiency indices 



Fig. 8.1. Errors, error estimates, and summed error indicators (top) and efficiency indices 
(bottom) versus the number of degrees of freedom N for adaptively refined meshes for the Poisson 
problem. Note the excellent agreement between the error rj (dashed black curve) and the error 
estimate rjh (solid red curve), as well as the convergence of the efficiency index rjh/rj towards 1. 



We consider the discretization of these equations by a mixed finite element space 
Vh = Vh consisting of the tensor fields composed of two first-order Brezzi-Douglas- 
Marini elements for the stress tensor, piecewise constant vector fields for the displace- 
ment, and continuous piecewise linears for the rotation [151 [16]. 

We consider the domain ft = (0,1) x (0,1) and the exact solution u{x,y) = 
{xy sm{7ry) ^ 0) for /i = 1 and A = 100, and insert 

._i . X / 7tll(2x cos(7ry) — 7Txysm(7Ty)) 

^ ^ y /j.[7ry cos[7ry) -\- sm[7ry)) -\- A[7ry cos[7ry) -\- sm[7ry)^ 

As a goal functional, we take a weighted measure of the average shear stress on the 
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right boundary, 

M((cr,ii,7)) = J cr-n-(V^,0) d5 ^ -0.06029761071, 

where F = {(x, y) : x = 1} and ip = y{y — 1). 

The resulting errors, error estimates, error indicators, and efficiency indices are 



plotted in Figure |8.2[ Again, we note that the error estimate rjh is very close to 
the actual error r]. We also note the good performance of the error indicators that 
overestimate the error by around a factor of 2 — 4. This is remarkable, considering that 
that the error estimate and error indicators are derived automatically for a non-trivial 
mixed formulation and involve automatic extrapolation of the dual solution from a 
mixed [BDMi]^ x DGq x Pi space to a mixed [BDM2]^ x DGi x P2 space. As far as 
the authors are aware, this is the first demonstration of goal-oriented error control for 



this discretization of the formulation (8.5) 



8.3. The stationary incompressible Navier— Stokes equations. Finally, 
we consider a stationary pressure-driven Navier-Stokes flow in a two-dimensional 
channel with an obstacle. We let Q = Qc\^o^ where Qc = (O?^) x (0,1) and 
no = (1.4, 1.6) X (0, 0.5). We let = {{x, y) e dn,x = or x = 4} denote the Neu- 
mann (inflow/outflow) boundary and let = 1] \ Qn denote the Dirichlet (no-slip) 
boundary. 

We consider the following nonlinear variational problem for the solution of the 
stationary incompressible Navier-Stokes equations: find (u^p) G V such that 

F{{u,py,{v,q)) = 

for all {v^q) G where 

F{{u, p); {v, q)) = i^(grad u, grad v) + (grad u-u.v)- {p, div v) + (div u, q) + (pn, v)dnN • 

Here, p is a given boundary condition at the inflow/outflow boundary. 

The trial and test spaces are given hy V = V = Hi q^^{VL]M?') x L'^{Q). We let 
the (kinematic) viscosity he v = 0.02 and take p = 1 at x = (inflow) and p = at 
X = 4 (outflow). The quantity of interest is the outflux at x = 4, 

M(u,p)= i ii-nds?^ 0.40863917. 

Jx=A 

The system is discretized using a Taylor-Hood elements; that is, the velocity space is 
discretized using continuous piecewise quadratic vector fields and the pressure space 
is discretized using continuous piecewise linears. The nonlinear system is solved using 
a standard Newton iteration. 



The results for this case are shown in Figure 8.3 As seen in this figure, the error 
estimate is not as accurate as for the two previous test cases. The efficiency index os- 
cillates in the range 0.2 — 1.0. This is not surprising, considering that (i) a linearization 
error is introduced when linearizing the dual problem around the computed approxi- 



mate solution rather than computing the average (6.3), and (ii) both the primal 
and dual problems exhibit singularities at the reentrant corners making the higher- 
order extrapolation procedure suboptimal for approximating the exact dual solution. 
Still, we obtain reasonably good error estimates and error indicators. Furthermore, 
the adaptive algorithm performs very well when comparing the convergence obtained 
with the adaptively refined sequence of meshes to that of uniform refinement, cf. Fig- 
ure 8.4 The final mesh is shown in Figure 8.5, and we note that it is heavily refined 



in the vicinity of the reentrant corners. 
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Fig. 8.2. Errors, error estimates, and summed error indicators (top) and efficiency indices 
(bottom) versus the number of degrees of freedom N for adaptively refined meshes for the mixed 
elasticity problem. Note the excellent agreement between the error rj (dashed black curve) and the 
error estimate r]h (solid red curve), as well as the convergence of the efficiency index ijh/v towards 1. 



9. Conclusions. We have demonstrated a new strategy for automated, adaptive 
solution of finite element variational problems. The strategy is implemented and freely 
available as part of the DOLFIN finite element library [29]; accessible both through 
the Python and C++ interfaces. 

The strategy and its implementation are currently limited to stationary nonlinear 
variational problems. Another limitation is the restriction to conforming finite element 
discretizations. These are both issues that we plan to consider in future extensions 
of this work. Additionally, we have assumed that the dual problem is well-posed. 
This assumption may fail in cases where the primal problem has been stabilized by 
the introduction of additional terms; the adjoint (linearized) dual problem is then not 
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(a) Errors 




(b) Efficiency indices 



Fig. 8.3. Errors, error estimates, and summed error indicators (top) and efficiency indices (bot- 
tom) versus the number of degrees of freedom N for adaptively refined meshes for the Navier-Stokes 
problem. This is a detail of the full convergence plot shown in Figure where the convergence of 
the adaptive algorithm is contrasted to the convergence obtained with uniform refinement. 



necessarily well-posed. Automated error control for such formulations is an interesting 
topic for further research, but is beyond the scope of the present work. 

Although the implementation has been tested on a number of model problems 
with convincing results, the effect of the linearization error (approximating u ^ Uh 
in (6.3)) is unknown. As a consequence, the computed error estimates typically un- 
derestimate the error for nonlinear problems. The effect of the linearization error and 
its proper treatment remains an open (and fundamental) question. Also, the extrap- 
olation algorithm proposed and numerically tested here should be examined from a 
theoretical viewpoint. 

We remark that the techniques described in this paper could also be used for 
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norm-based error estimation. A posteriori error estimates for energy or other Sobolev 
norms typically rely on computing appropriately weighted norms of cell and averaged 
facet residuals. Hence, the strategy described here provides a starting-point for the 
automatic generation of norm-based error estimators. 
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